Sequential design with early stopping (restricted action set) - risk based
Published
May 19, 2025
Modified
May 14, 2025
Load simulation results
# Each input file corresponds to the results from a single simulation# scenario/configuration.# Load all the files into a single list.# files of interestsim_lab <-"sim05-15"# files of interestflist <-list.files(paste0("data/", sim_lab), pattern ="sim05")toks <-list()l <-list()i <-1for(i in1:length(flist)){ l[[i]] <- qs::qread(file.path(paste0("data/", sim_lab), flist[i])) toks[[i]] <-unlist(tstrsplit(flist[i], "[-.]"))}ix_scenario <-9ix_trial =599
Results from example trial
Cumulative probability of each decision type
# Cumulative probability of decisions:# Traverse the list of simulation results and for each one summarise the # cumulative probability of each decision type.i <-1d_tbl_1 <-data.table()# For each scenario that was simulatedfor(i in1:length(l)){# extract the decision matrix - sim, analysis, quantity, domain level decision d_dec_1 <-copy(l[[i]]$d_decision)# config for scenario l_cfg <-copy(l[[i]]$cfg)# number of enrolments at each interim (interim sample size sequence) d_N <-data.table(analys =seq_along(l_cfg$N_pt), N = l_cfg$N_pt)# long version of decisions d_dec_2 <-melt(d_dec_1, id.vars =c("sim", "analys", "quant"), variable.name ="domain")# Should be right, but just in case...if(any(is.na(d_dec_2$value))){message("Some of the decision values are NA in index ", i, " file ", flist[i]) d_dec_2[is.na(value), value :=FALSE] } d_dec_2[, domain :=as.numeric(gsub("d", "", domain))]# Domains 1, 3 and 4 will stop for superiority or futility for superiority.# Domaain 2 will stop for NI or futility for NI.# No other stopping rules apply and so we only evaluate the operating # characteristics on these, i.e. we do not care about the results for the # cumualative probability of ni for domain 1, 3 and 4 because we would never# stop for this. d_dec_2 <-rbind( d_dec_2[domain %in%c(1, 3, 4) & quant %in%c("sup", "fut_sup")], d_dec_2[domain %in%c(2) & quant %in%c("ni", "fut_ni")] )# compute the cumulative instances of a decision being made by sim, each # decision type and by parameter d_dec_2[, value :=as.logical(cumsum(value)>0), keyby = .(sim, quant, domain)] d_dec_2 <-merge(d_dec_2, d_N, by ="analys")# cumulative proportion for which each decision quantity has been met by # analysis and domain d_dec_cumprob <- d_dec_2[, .(pr_val =mean(value)), keyby = .(analys, N, quant, domain)] d_tbl_1 <-rbind( d_tbl_1,cbind(scenario = i, desc = l_cfg$desc, d_dec_cumprob) )}
Expected number of enrolments
# Similar to above but focus on expected number of enrolments# Traverse the list of simulation results and for each one summarise the # sample sizes at which stopping for a domain occurs for any reason.# All we are trying to get to is the expected sample size by domain and # are not really interested in what decision was made. The cumulative prob# of each decision type is computed previously.i <-1d_tbl_2 <-data.table()for(i in1:length(l)){# extract the decision matrix that contains the stopping decision by # quantity (superiority, ni, futility, etc) for each domain by analysis d_dec_1 <-copy(l[[i]]$d_decision)# config for scenario l_cfg <-copy(l[[i]]$cfg)# interim looks d_N <-data.table(analys =seq_along(l_cfg$N_pt), N = l_cfg$N_pt)# long version of decision, e.g.# sim analys quant domain value# <int> <int> <char> <fctr> <lgcl># 1: 1 1 sup d1 FALSE# 2: 1 2 sup d1 FALSE# 3: 1 3 sup d1 FALSE# 4: 1 4 sup d1 FALSE# 5: 1 5 sup d1 FALSE# 6: 2 1 sup d1 FALSE d_dec_2 <-melt(d_dec_1, id.vars =c("sim", "analys", "quant"), variable.name ="domain")# Should be right, but just in case we had a sim fall over, fill in valueif(any(is.na(d_dec_2$value))){message("Some of the decision values are NA in index ", i, " file ", flist[i]) d_dec_2[is.na(value), value :=FALSE] } d_dec_2[, domain :=as.numeric(gsub("d", "", domain))] d_dec_2 <-merge(d_dec_2, d_N, by ="analys")# Extract the first instance of any decision occurring by sim and domain.# Domains 1, 3 and 4 will stop for superiority or futility for superiority.# Domaain 2 will stop for NI or futility for NI.# Sometimes a decision rule will not be hit for a domain and it will continue# to the max sample size. We will deal with that in a minute. d_dec_stop <-rbind( d_dec_2[ domain %in%c(1, 3, 4) & value == T & quant %in%c("sup", "fut_sup"), .SD[1], keyby = .(sim, domain)], d_dec_2[ domain %in%c(2) & value == T & quant %in%c("ni", "fut_ni"), .SD[1], keyby = .(sim, domain)] )setnames(d_dec_stop, "N", "N_stopped")setnames(d_dec_stop, "value", "stopped_early")# Add in any rows for which no early stopping happened d_dec_stop <-merge( d_dec_stop[, .SD, .SDcols =!c("analys")], # all combinations of sim and domainunique(d_dec_2[, .(sim, domain)]),by =c("sim", "domain"), all.y = T)# If domain or trial not stopped then record as having run to the # maximum sample size with no decision made. d_dec_stop[is.na(N_stopped), N_stopped :=max(d_N$N)] d_dec_stop[is.na(stopped_early), stopped_early := F]# So now we know where every domain was stopped and the reason for stopping# within each sim. Great. d_dec_stop[is.na(quant), quant :="null"] d_tbl_2 <-rbind( d_tbl_2,cbind(scenario = i, desc = l_cfg$desc, d_dec_stop ) )}
Distributions of posterior means (unconditional)
# Distribution of posterior means for parameters of interest.# Some simulated trials will have stopped prior to the maximum sample size and# these will have NA for their posterior means. If you were to summarise the # posterior means, they would thus be conditional on the trial having 'survived' # until the relevant interim. This means that you have missing data at later # interims, which creates a selection bias in that your selection of sims at any# given interim are not a random sample, but rather a sample conditioned on the # stopping rules. # If you do not account for this in some way then a summary can be either # optimistic or pessimistic depending on how the stopping rules interact # with the data. Here we try to account for this missingness by imputing the # missing posterior means with locf within each simulation.# Note that this is really only a partial 'fix' to get a sense of whether # what we estimate is representative of the parameter values we used to simulate# the data.i <-2d_tbl_3 <-data.table()for(i in1:length(l)){# config for scenario l_cfg <-copy(l[[i]]$cfg)# params d_pars <-copy(l[[i]]$d_post_smry_1) d_pars <- d_pars[par %in%c("lor", "rd")]# interim looks d_N <-data.table(analys =seq_along(l_cfg$N_pt), N = l_cfg$N_pt) d_pars <-dcast(d_pars, sim + id_analys + domain ~ par, value.var =c("mu", "se"))# locf d_pars[, `:=`(mu_lor =nafill(mu_lor, type ="locf"),mu_rd =nafill(mu_rd, type ="locf"),se_lor =nafill(se_lor, type ="locf"),se_rd =nafill(se_rd, type ="locf") ), keyby = .(sim, domain)]setnames(d_pars, "id_analys", "analys")# d_pars <-merge(d_pars, d_N, by ="analys") d_tbl_3 <-rbind( d_tbl_3,cbind(scenario = i, desc = l_cfg$desc, d_pars[, .(analys, sim, domain, N, mu_lor, mu_rd, se_lor, se_rd)] ) )}
Cumulative probability of each decision type
# > names(l[[ix_scenario]])# [1] "cfg" "d_pr_sup" "d_pr_ni" "d_pr_sup_fut" "d_pr_ni_fut" "d_decision" # [7] "d_post_smry_1" "d_all" "d_n_units" "d_n_assign" # extract the decision matrix - sim, analysis, quantity, domain level decisiond_N <-data.table(analys =0:5,N =seq(0, 2500, by =500))# Extract example trial specific datal_cfg <-copy(l[[ix_scenario]]$cfg)d_dat <-copy(l[[ix_scenario]]$d_all[sim == ix_trial])# Parameter summaries from analysisd_mod <-copy(l[[ix_scenario]]$d_post_smry_1[sim == ix_trial])# Decision matrix for this triald_dec <-copy(l[[ix_scenario]]$d_decision[sim == ix_trial])# Probability level reached by domain and analysisd_pr_sup <-copy(l[[ix_scenario]]$d_pr_sup[sim == ix_trial])d_pr_sup_fut <-copy(l[[ix_scenario]]$d_pr_sup_fut[sim == ix_trial])d_pr_ni <-copy(l[[ix_scenario]]$d_pr_ni[sim == ix_trial])d_pr_ni_fut <-copy(l[[ix_scenario]]$d_pr_ni_fut[sim == ix_trial])# observed datasetnames(d_dat, "id_analys", "analys")setnames(d_mod, "id_analys", "analys")d_mod <-merge(d_mod, d_N, by ="analys")d_pr_sup <-merge(d_pr_sup, d_N, by ="analys")d_pr_sup <-melt(d_pr_sup[, .(analys, d1, d3, d4, N)], id.vars =c("analys", "N"), variable.name ="domain")d_pr_sup[, domain :=as.numeric(gsub("d", "", domain))]d_pr_sup[, quant :="sup"]d_pr_sup_fut <-merge(d_pr_sup_fut, d_N, by ="analys")d_pr_sup_fut <-melt(d_pr_sup_fut[, .(analys, d1, d3, d4, N)], id.vars =c("analys", "N"), variable.name ="domain")d_pr_sup_fut[, domain :=as.numeric(gsub("d", "", domain))]d_pr_sup_fut[, quant :="sup_fut"]d_pr_ni <-merge(d_pr_ni, d_N, by ="analys")d_pr_ni <-melt(d_pr_ni[, .(analys, d2, N)], id.vars =c("analys", "N"), variable.name ="domain")d_pr_ni[, domain :=as.numeric(gsub("d", "", domain))]d_pr_ni[, quant :="ni"]d_pr_ni_fut <-merge(d_pr_ni_fut, d_N, by ="analys")d_pr_ni_fut <-melt(d_pr_ni_fut[, .(analys, d2, N)], id.vars =c("analys", "N"), variable.name ="domain")d_pr_ni_fut[, domain :=as.numeric(gsub("d", "", domain))]d_pr_ni_fut[, quant :="ni_fut"]d_pr_dec <-rbind(d_pr_sup, d_pr_sup_fut, d_pr_ni, d_pr_ni_fut )d_pr_dec[domain %in%c(1, 3, 4), question :="Superiority"]d_pr_dec[domain %in%c(2), question :="Non-inferiority"]# domain 1 relates only to the late acute silod_dat_d1 <- d_dat[ silo ==2, .(y =sum(y), n =sum(N)), keyby = .(silo, analys, d1)]d_dat_d1[, `:=`(cy =cumsum(y), cn =cumsum(n)), keyby = .(silo, d1)]d_dat_d1[, p_obs := cy/cn]d_dat_d1[, d1 :=factor(d1, labels =c("DAIR", "One-stage", "Two-stage"))]# domain 2 relates only to the group that receive one-staged_dat_d2 <- d_dat[ d1 ==2, .(y =sum(y), n =sum(N)), keyby = .(silo, analys, d2)]d_dat_d2[, `:=`(cy =cumsum(y), cn =cumsum(n)), keyby = .(silo, d2)]d_dat_d2[, p_obs := cy/cn]d_dat_d2[, d2 :=factor(d2, labels =c("Selected", "12wk", "6wk"))]# domain 3 relates only to the group that receive two-staged_dat_d3 <- d_dat[ d1 ==3, .(y =sum(y), n =sum(N)), keyby = .(silo, analys, d3)]d_dat_d3[, `:=`(cy =cumsum(y), cn =cumsum(n)), keyby = .(silo, d3)]d_dat_d3[, p_obs := cy/cn]d_dat_d3[, d3 :=factor(d3, labels =c("Selected", "none", "12wk"))]# domain 4 relates only to the group that receive two-staged_dat_d4 <- d_dat[ , .(y =sum(y), n =sum(N)), keyby = .(silo, analys, d4)]d_dat_d4[, `:=`(cy =cumsum(y), cn =cumsum(n)), keyby = .(silo, d4)]d_dat_d4[, p_obs := cy/cn]d_dat_d4[, d4 :=factor(d4, labels =c("Selected", "no-rif", "rif"))]d_dec <-melt(d_dec, measure.vars =paste0("d", 1:4), variable.name ="domain")d_dec[, domain :=as.numeric(gsub("d", "", domain))]# Subset to the decisions that are evaluated by domaind_dec <-rbind( d_dec[domain %in%c(1, 3, 4) & quant %in%c("sup", "fut_sup")], d_dec[domain %in%c(2) & quant %in%c("ni", "fut_ni")])d_dec_thres <-data.table(quant =rep(c("sup", "sup_fut", "ni", "ni_fut"), each =4),domain =rep(1:4, len =16) )d_dec_thres[, threshold :=c( l_cfg$dec_probs$thresh_sup, l_cfg$dec_probs$thresh_fut_sup, l_cfg$dec_probs$thresh_ni, l_cfg$dec_probs$thresh_fut_ni )]## State transition through the different states of knowledge based on decisionsd_dec_timeline <-merge(CJ(domain =1:4, analys =0:5), d_N, by ="analys", all.y = T)setorder(d_dec_timeline, domain, analys)d_dec_timeline[domain %in%c(1, 3, 4), question :="Superiority"]d_dec_timeline[domain %in%c(2), question :="Non-inferiority"]d_dec_timeline[N ==0, decision :="Indeterminate"]i <- j <-1# For i across domainsfor(i in1:4){# For j across analysesfor(j in1:5){if(i %in%c(1, 3, 4)){if(d_dec[domain == i & analys == j & quant =="sup", value]) { d_dec_timeline[domain == i & analys == j, decision :="Superior"] } elseif (d_dec[domain == i & analys == j & quant =="fut_sup", value]) { d_dec_timeline[domain == i & analys == j, decision :="Futile (sup)"] } else { d_dec_timeline[domain == i & analys == j, decision :="Indeterminate"] } } else {if(d_dec[domain == i & analys == j & quant =="ni", value]) { d_dec_timeline[domain == i & analys == j, decision :="NI"] } elseif (d_dec[domain == i & analys == j & quant =="fut_ni", value]) { d_dec_timeline[domain == i & analys == j, decision :="Futile (NI)"] } else { d_dec_timeline[domain == i & analys == j, decision :="Indeterminate"] } } }}d_dec_timeline[, decision :=factor(decision, levels =c("Indeterminate", "Superior", "Futile (sup)", "NI", "Futile (NI)"))]# Dealing with those for which a decision was not reached.d_decision <-data.table(domain =1:4,N =2500)if(any(d_dec_timeline$decision !="Indeterminate")){ d_decision <- d_dec_timeline[ decision !="Indeterminate", .(N =min(N)), keyby = .(domain, decision)]}d_decision <-merge( d_decision, data.table(domain =1:4),by ="domain", all.y = T)d_decision[is.na(decision), `:=`(decision ="Intermediate", N =2500)]d_dec_timeline <- d_dec_timeline[N <=max(d_decision$N)]
Table 1 shows the decisions made for each domain (or indeterminate if no decisions were made).
Figure 1 shows the knowledge transitions based on the decisions made for each domain by sample size up to the point where the trial was stopped either due to running out of resources or having addressed all the questions.
Initially, all domains start in an indeterminate state in that neither treatment arm is preferred. As the data accrues and analyses progresses, the knowledge state for each domain may transition to, superiority, non-inferiority or futility.
Code
p1 <-ggplot(d_dec_timeline, aes(x = N, y = decision)) +geom_point() +scale_y_discrete("", drop=FALSE) +facet_wrap(domain ~ question, labeller = label_both)suppressWarnings(print(p1))
Of those receiving one-stage (see figure for domain 1) and irrespective of silo membership, approximately 70% are assumed to enter the AB duration domain.
Of those receiving two-stage (see figure for domain 1) and irrespective of silo membership, approximately 90% are assumed to enter the AB ext-proph domain.
Across the entire trial sample, approximately 60% are assumed to enter the AB choice domain.
Here, we know there is an effect of AB choice and would hope that a decision for superiority is made such that new participants are directed to receive rifampacin.
Figure 8 shows the probability associated with each decision type for the randomised comparisons by domain and enrolment progression.
For example, for domain 1 in analysis 1, the probability of revision being superior to DAIR (per the definition of superiority, which is a measure that is relative to nominated reference level for superiority) is approxiamtely 0.55. To make a superiority decision, this probability needs to exceed 0.920.
Similarly, for domain 1 in analysis 1, the probability of the superiority decision being futile (per the relevant definition) is approximately 0.30. To make a futility decision (related to superiority) this probability has to fall below 0.300.
Note:
The futility probabilities are based on being below a given threshold (rather than above).
Code
p1 <-ggplot(d_fig_1, aes(x = N, y = value)) +geom_point(aes(col=quant), position =position_dodge2(width =100), size =0.6) +geom_linerange(aes(ymin =0, ymax = value, col=quant), position =position_dodge2(width =100), lwd =0.25)+geom_hline(data = d_fig_2,aes(yintercept = threshold, col=quant), lwd =0.25 ) +scale_x_continuous("") +scale_y_continuous("Decision probability", breaks =seq(0, 1, by =0.2)) +scale_color_discrete("") +facet_wrap(domain ~ question, labeller = label_both)suppressWarnings(print(p1))
Figure 8: Probability decision summaries
Reveal scenario
The results are from a simulated trial picked from scenario 9: Null effect in all domains +silo specific effects.